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Abstract 

In this paper we present the numerical analysis of spectral methods when non-constant 
coefficients appear in the equation, either due to the original statement of the equations 
or to take into account the deformed geometry. A particular attention is devoted to the 
optimality of the discretization even for low values of the discretization parameter. The 
effect of some “overintegration” is also addressed, in order to possibly improve the accuracy 
of the discretization. 
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1. INTRODUCTION. 


Spectral methods were introduced around 20 years ago by S. A. ORSZAG in order to provide a high 
order accuracy for the numerical simulation of partial differential equations. The original idea was to use 
truncated Fourier aeries to approximate the (smooth) solution when the problem was supposed to be provided 
with periodic boundary conditions. In order to tackel problems with more general boundary conditions 
(Dirichlet or Neumann type), the set of (algebraic) polynomials replaced the set of truncated series, but the 
characterization of the unique discrete function that would provide the numerical solution was still achieved 
following the original strategy. This enters in the context of the collocation method where the numerical 
solution is chosen so as to satisfy the original partial differential equation at some suitably chosen collocation 
points. Of course the choice of the set of collocation points is of fundamental importance for the accuracy 
of the method and the first remark is to notice that the number of collocation points must be equal to the 
dimension of the space of approximation. Otherwise, the problem could, in general, be overspecified. 

When other authors arrived to the analysis of the problem, ([1], [2], [3], [4] among others) they realized 
that the collocation method could be interpreted as a variational problem with numerical integration. In 
fact, it is quite common now to realize that the spectral methods are very close to the finite element method 
in its p or h — p version, where convergence is achieved by increasing the order of the polynomial degree and 
not by diminishing the size of the elements. It is in this framework that the domain of application of the 
spectral methods has been generalized. The plain spectral method suffers from being constrained to very 
simple domains; they are limited to be slightly deformed squares (in 2-D) or cubes (in 3-D). The idea to 
couple domain decomposition techniques to the spectral discretization developed rapidly in order to cope 
with this initial drawback ([5], [6]). However, by starting from the strong formulation of the equations, this 
can only produce the Schwarz algorithm ([6], [7]) or a strong coupling between the elements ([5]), where the 
solution (in the case of second order elliptic problems) is searched as a global C l function that is piecewise 
polynomial. This results in some drawbacks. First, in the Schwarz algorithm, there is an increase of the 
work due to the double computation over the overlapping region (recall that this one has to be large enough 
in order to achieve a good convergence rate of the algorithm [7], [8]). Second there is a lack of optimality of 
the approximation in the strong coupling formulation — both from the numerical analysis point of view and 
from the algorithmic point of view — a consequence of an overconstraint problem due to the C l matching. 

In this context, understanding of the similarity between the collocation method and a variational for- 
mulation used with consistent numerical quadrature brings a lot of flexibility. Indeed, the coupling required 
in the nonoverlapping decomposition of the domain is weaker (only C°) and allows for constructing an op- 
timal method. This hits lead to the spectral element method [9], [10], [11] and more recently to the mortar 
element method [12], [13], [14]. The variational method involves integrals that can be computed (with or 
without numerical quadratures) separately over each subdomain. The spectral element method uses consis- 
tent quadrature formulas and in this respect (but not only this one) conserves the spirit of former spectral 
discretization (the other points being the use of tensorial basis and tensorial evaluations of the residuals). It 
has been shown ([10], [15]) that the method can still be interpreted as a collocation method within each ele- 
ment wheras a suitable equation is satisfied at the interface of each subdomain. There are still fundamental 
differences between the spectral and the finite element methods that are important as regards the numerical 
implementation of the method, however the general philosophy is the same. The point that we want to 
address in this paper is related to one of the differences between the p or h — p version of the finite element 
method and the spectral methods, more precisely to one possible drawback of the spectral method. Indeed, 
as we have said, the spectral method can be interpreted in a finite element framework when a particular nu- 
merical integration formula is used. Derived from the collocation method, the numerical integration formula 
is constrained to be based on a certain number of points. This number is not related to the fact that some 
quantity must be well (or exactely) computed, but has to be equal to the dimension of the discrete space. 
Using a vocabulary that is standard in the finite element context, the points of the numerical quadrature 
formula have to be unisolvent with respect to the discrete space. This is a restriction with respect to the fact 
that in finite element methods, people use more general integration formulas based on accuracy considera- 
tions. In this sense, there is often an overintegration in the finite element statement. The current spectral 
methods use these “consistent” (with respect to the numerical discretization) quadrature rules with a priori 
no lack of accuracy. As shown in the appendix, the possibility to use an overintegration in the context of 


1 


spectral methods is however possible, but is not of importance in most interesting cases as demonstrated in 
this paper. 

The situations where overintegration could be of importance are those where non-constant coefficients 
are involved, either because they are present in the original formulation of the equations or because they are 
introduced when a deformed geometry is considered and a mapping to a domain of reference is used. Indeed, 
it is in these cases that the numerical integration based on N + 1 Gauss-Lobatto points is not accurate enough. 
We analyze in details these two situations and point out the cases where overintegration is required. The 
analysis is illustrated with many numerical experiments. This allows to strengthen our theoretical analysis 
since the exact value of some constants could be of importance and are not exactely taken into account in 
the theoretical analysis. 

The analysis of the effect of overintegration has never been addressed in the spectral context. Previous 
analysis of Legendre spectral approximations of problems with non-constant coefficients [10] or with deformed 
geometries [7] could in this sense be misleeding as the results presented there were proving a spectral type 
convergence (i.e. faster than any algebraic rate if the solution is analytic) but were not optimal. The basic 
ingredient to get optimality is a new result of [16] and deals with the fact that the interpolation at Gauss-type 
points is optimal. 

The paper is organized as follows: 

In section 2, we consider the case of original non-constant coefficients equations and we analyze in which 
cases the overintegration may improve the accuracy of the approximation to allow for optimal results. An 
appendix presents some numerical considerations on the implementation that shows that the tensor product 
evaluations can still be preserved even in the case of overintegrations. 

In section 3, we analyze the case where the non-constant coefficients are induced by the treatment of a de- 
formed geometry. We prove that in this situation the overintegration is unnecessary to provide the optimality 
of the approximation. 

In sections 2 and 3 we also present some remarks related to the impact of our analysis as regard the design 
of the best schemes in the case of deformed geometries. 


Acknowledgements: The question of the importance of overintegration has been raised during a meeting 
with I. BABUSKA. The many discussions with A. T. PATERA have been very enlightening. We thank both 
of them for their help. This research has been supported in part by NSF under ASC-8806925, ONR and 
DARPA under contract N00014-85-K 0208. Part of this work has been done while the first author was in 
residence at ICASE, NASA Langley Research Center, Hampton, VA. 
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2. ANALYSIS OF SPECTRAL APPROXIMATION WITH NON CONSTANT COEFFICIENTS. 

2.1 — The One Dimensional Case. 

Let us consider the following problem: Find u 6 #o(— 1, 1)* such th&t 

-(p«.). + }« = /, in (-1,1). (1) 

This problem is supposed to be well posed in Hq( — 1, 1) in the sense that the bilinear form a, defined over 
(Ifo(— 1, l)) 3 as follows 

a(u,v) = J pu m v z dx + J quvdx y (2) 

is elliptic and continuous over ( jETo ( — 1, l)) 3 . In fact, for sake of simplicity, we shall assume that there exists 
three constants Pi,Pa and q 2 such that 

0 < Pi < P(») < Pa. Vx 6 (-1, 1); (3) 

0 <$(*)< $2. Vx G (—1, 1). (4) 

In addition, the forcing function / is supposed to be at least in L 2 (— 1, 1). We want to approximate this 
problem by a spectral element method. Our first work will be to design the discrete space that will approx- 
imate J/q(— 1, 1). Given a fixed integer K } we shall first consider a partition of ( — 1, 1) in K subintervals /*, 
where 

Ik = (ojk,a 4+ i), (5) 

and the a* are ( K -f 1) points in (-1,1) such that 


— 1 = oq < «i < . . . < a/f-i < a* = 1, 


(and the lengh ot+i — a* is supposed to be 0(K *).) 

The spaces of discretization are imbedded in 

y n = {* el 2 (-l,l)i € 2P*(/*)}, (6) 

where N is some integer and JP^(A) denotes the set of all polynomials of degree < N over A. This space 
will discretize 2/ 3 (— 1, 1). In order to discretize Hq(— 1, 1), let us introduce 

X N = Y N n HU-1,1), ( 7 ) 

(which means that ^|/ fc (a*+i) = ^|/ fc+l (ok+i) and $ vanishes at ±1.) 

The second step in the discretization process iB to define the discrete problem. As indicated in the 
introduction, the point of departure is the following continuous variational formulation of the problem (1): 
Find u G Af<J(— *» 1) BUC ^ that 

a(u,v) = (/,v), Vv e !). ( 8 ) 

where the notation (.,.) stands for the L 3 (-l, l)-scalar product. Given an integer Af, M > N , the discrete 
problem is: Find u# G such that 

a>M (uat, vn) = (/,v^)m , E Jfjv* (9) 

The definition of clm and (.,.)*/ from a and (.,.) will use a composed Gauss-Lobatto quadrature formula. 
More precisely, this one is based on the data, over each segment J* of a set of M + 1 points f * 

= (o < < ••• < £m-1 < iu ~ °k+l« ( 10 ) 
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and a set of associated weights p* that are such that 


KM * 

5> S £5>(tf)p? = / *(*)<**, V4> 6 Y 2 U- 1 . ( 11 ) 

GL *=0 i=0 Jl 

The definition of au and i* now as follows: 

au(^N,v N ) - Y^PVNzVn* +'%2quN'>N, V(u N ,v N ) e Xj,, (12) 

GL GL 

(*,*)* = v(^)ey| (13) 

Remark 2.1. The standard spectral element method utilizes a consistent Gauss-Lobatto quadrature rule, 
i.e. M = N. The case M > N corresponds to an overintegration with respect to the standard spectral 
element method. Such an overintegration might be of interest to introduce in the case of non-constant 
coefficients. The analysis that follows indicates the situations for which this is the case. 

The interesting feature about (., .)m i» that it defines over Ym a discrete scalar product that is uniformly 
equivalent to the L 2 ( — 1, 1) one. Indeed, it has been proven in [4] that 

(*, <t>) < (4>, 4)u < 3 (*, *), V<t> 6 Y U . (14) 


The following Lemma is then completely standard 

Lemma 2.1. There exists one and only one solution ujv to problem (9) and there exists a constant C such 


that 


Proof: The Lemma follows easily from the constatation that um is uniformly continuous and elliptic over 
X# t i.e. there exists two positive constants a and 7 such that 


<o Af (u A r,UAf) <7|M|tf.(_x,i) (15) 

and can be deduced from (3), (4) and (14). 


Our further analysis is devoted to the derivation of optimal error bounds for the numerical solution. 
The following Lemma is also standard and follows from (15). 

Lemma 2.2. There exists a constant C independent of N such that 


' WNtXu v N iX„ ||vjV||tf>(-l,l) 

+ .up (/ ',T' ,) ,r </ ’ "")" ]• 


(16) 


The following Lemma proves that the last term in (16) is of the same order as the best L a (— 1, l)-fit of 
/ by polynomials of degree 2Af — N. Let us first denote by C^-(— 1, 1) the space 


C° K (- 1 , 1 ) = € L\- 1 , 1 ), 3 ^* € C°(f k ),^ Ih = rj> h } 

and similarly, for any n, by H^(- 1, 1) the space 

*£(-l.l) = {*€l a (-M), 4>\l„€H^I k )}, 


(17) 

(18) 


we have 


4 



Lemma 2.3. Let p be a real number > 1/2. For any $ £ 1), we have 


sup 

vu£X*r 


kw 6 v£-» ||vaHL’(- 1 , 1 ) 


i c < 2 " - *>- 1 


( 19 ) 


( 20 ) 


i=l 


Proof: Let us first consider the case where AT = N. Let us introduce, for any integer L, the L 2 (— 1, 1)- 
projection V’i = Il£,(^) of ^ over 2Px,(— 1, 1), we have 


(V'. t'w) - (iMa r)N = [(iMw) - (V’w.vw)) + [('Pn^n) - ('Pn,v n ) n ] 

+ [(‘•pNtV^N - ('PiVn)n]- 


(21) 


It » an easy matter to note that 

KV’.vjv) - (V’aTiVjv)I < C\\i> - V’iv|U>(-i > i)||vtflU , (-M)* ( 22 ) 

Besides, let us denote by 7jv the operator of interpolation from the set C^(— 1, 1) onto Yu such that 

= Vi = 0 N, Vk=l,...,K. (23) 


then 


KV’at, vn)n - U>,vn)n | = KV’/r -Zn'!>,'>n)n\ 

< (ipN - 1 n1>, *I>N - ^N^>)n 2 (vu, VN )J/ 3 , 


which, thanks to (14), gives 


KV’at. vat)at - (V** vat)at| < C\\il>N -7^V’IU’(-i,i)l|t'Ar|U>(-i 1 i) 

< ciW'P - V’wlU^-i.i) + IIV 1 - ^wV’IU»(-i,i))ll"wlU > (-i,i)- ( 24 ) 

It remains to estimate the middle term in (21). Let us first assume that K = 1, corresponding to a plain 
spectral method with no domain decomposition. The term in question can then be written as follows 

{'l>N> v n) — {'Pni vn)n — ($wvjv)[(-Lw> Ln ) — (Ln, Ln)n] 


where, for any j € IN 


(Li.LiY ’ (Lj.Lj) 1 


(25) 


are the coefficients of rp (and thus of i>n) and v N respectively in the basis of the Legendre polynomials 
It is then an easy matter to note, from (14), that 


KV’tf.VAr) - (V’jv, vn)nI < C|^//®w|||Ljv||!>(_j > i) 

< C\j> N \ ||ijv||£,>(-l,l)|t)jv| ||iAr|U>(-i,i) 


Also, it follows easily that 

HMHMUm-M) < 11^ - V’rv-i|U 3 (-i,i) 

and 

l®Jv||| jC 'Ar||i,a(_l,l) < ||VAT - HjV— l(vjsr )!|i a ( — 1 . 1 )> 

so that we deduce 

I(V’jv.iw) ~ (^jv.vjvJatI < C|| V* “ V’Ar-i||i>(-i,i)||t>jv - nAf_i(v w )|U»(-i,i)- ( 26 ) 
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The multidomain case where K > 1 is treated similarly over each subdomain J* and the error is the same 
(the notations are in this case just more complicated). Let us recall that we have for any p > 0, and any L 

II* - n*rf|U»(-t,i) < CL-» IMIh, ( _ m) , v* G 1 ). ( 27 ) 

we derive by summarizing (21), (22), (24), and (26) that 


(*■ vn) - 

vwe?*, ||«rv||tf»(-r,i) 


< 2||* - *w|U»(-i,i)+||* - 2//*IU J (-i,i) 

+ CN~ 1 \\ip - ^w-i||i-»(-i,i) 


This completes Lemma 2.3 in the case M - N, after we recall that for any p > 1/2, and any L, we have [16] 


II* - 2**IU»(-i,i) < C^ -P ll*lltf'(-M), V* € H»(- 1, 1). (28) 

so that the contribution of - V*Ar- 1 ||r,»{_ x,i) »« much smaller than the other terms. 

Let us now consider the case where M > N and let us use the same strategy as before, with 1 P 2 M-N in 
place of V’/Vi where V t > 2 M-n >® some good approximation of * of degree < 2 M — N that will not be necessarily 
its 2/ 3 (-l, 1) -projection and will be precised in the following. In order to prove Lemma 2.3, we must now 

estimate ... 

(lp,V N ) - (ip, V N ) M = (lp,V N ) - (lp 2 M-N,VN) + (*2M-JV,*»Ar) - (V >2M-n,vn)m 

+ (^>2M-N,Vn)u - (tP,Vn)m 

The first term and the second term are treated exactely in the same way as previously and we obtain 


KV’.t'tf) - (i> 2 M-N,VN)\ < II*- * 2 M-iv|U>(-l 1 l)||«Af||L>(-l,l)» 


( 30 ) 


and 

(^2M-N,Vn) - (-ip2M-N,t)N)M = iJ>2M-NVn[(L2M-N , L N ) - (L 2 M -N , L N ) M ]- 

£From the basic properties of the Legendre polynomials we recall that that 


(L m ,L ■„) = 


2^m,n 

2n + 1 ’ 


V(m,n), 


and that 

(n + l)L n+1 (x) = (2 n + l)xL n (i) - nl n _ i(x). 
A reiterate use of these two formulas yields 


/r r , 4M — 2N — 1 N+ 1 /r r , _ 

(L2M-n,Ln)m = — TU~~Tt — 2N + 1 ( l ™-n -i, ^>n+i)m 


2M -N 

< (Lm,Lm)m- 


(31) 


Now using (14), we derive that 

(LaM-Afi Ln)m < 3 I|7'm||l>(-i,i) ^ 3 I|7'2M-at||i,»(-i,i)||7'Jv|U»(-i,i). 

so that, as before, 


\(lp2M-N,VN) - (0>2M-N,Vn)m\ < CN X ||* - *2M-W-l|U*(-l,l). (32) 

The treatment of the last term in (29) has to be different. Indeed the same proceedure as before would result 
in a bound of this term by the quantity ||* — 1u V'||z, a (— x,i) an ^ thus a loss with respect to the optimality we 
could expect from the other terms. Besides, the numerical experiments (see fig 2.1) clearly shows a behaviour 
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like || ip — V'aAf-ArlU^-i.i)- In order to prove this type of convergence, let us go back to the case K = 1 and 
note that (28) yields, for any n > 1/2 

||Zj#*||l.(-i,i) < IMU’(-m) + M-'IMUh-m), V* € 1). (33) 

i,From (14) it now follows that 

( lf>3M-N , Vn)m - {'P, Vn)m = (V^M-AT - V>, V N ) M 

< 3||Iw(^ajif-w -^)||4»(-i,i)ll®Ar|U»(-i,i) 

< 3[||^aAf-AT - V’IU»(-1,1) + m-n - V’IIh<‘(-i,i)]|I w wIU , (-».i)- 

At this point the correct rp 2 M -n will be chosen. We see that it has to be a good approximation of xp in the 
2P(-1, l)-norm as well. This will be the case of the projection of rp with respect to the H*{- 1, l)-norm 
over Ytf. It is proven in [17] that this element satisfies 

||V» - ^2JM-jv||hj(-i,i) < Vi ' < n, (34) 

and from (29), (30) and (32), the proof of Lemma 2.3 is complete also in the case where U > N. 

The second inequality (20) is quite trivial since the middle terms in (21) and (29) simply vanish. 


We are now in position to Btate the following 

Theorem 2.4. Let us suppose that the solution u to problem (1) belongs to ifjf(— 1, 1), that p belongs 
to Hk(- 1, 1), q belongs to #£(-1, 1) a nd that f belongs to //* (-1, 1) and in addition that the four real 
numbers a - p, v and p are larger than 1/2. Then the following error estimate holds 


lit* - u*||*.(- M ) <C7[IV 1 '" <, ||u||tf»(_i i i) + (2 M - JV )-m<»(^..-i.rt [ || p ,| HS( _ 1(l) || tt ||^ +t ,, )( _ iii) 

+ + ll/lltf£(-l,l)]]- 

Proof: It is a direct application of Lemma 2.3 to derive 

8Up (/’ tw) - (/.^k < (2 M - Ar)-'’||/|| < ( _ M) . 

v N ex N ||wAf||tf‘(-i,i) 

Let us choose wn as being the best fit of u in the hf l (-l, l)-norm. It follows from [17] that 

||« - uw||tf‘(-i,i) < C7JV l_ "||t*||frj(_i,i). 

If we now let v>n be any polynomial in Xs , it follows from (3) and (4) that, 


(35) 


(36) 


(37) 


\a M (w N , V N ) - a M (w N , v w )| < C(p 2 , 92 )||t"JV - u>w||H>(-i,i)||i>/v||i/*(-i 1 i) 

< C(P2.92)[||w - tsjv||//>(-M) + ll u_ *Ar||i/»(-i 1 j)]||vw||i/‘(-».i) 


o(u,VAr) - OM(^Af>t>jy) < o(lI>Af,Vjv) -QSf(u>^,V/v) 

.JeL ||vw||i/‘(-i,i) “ ||vAr||tf»(-i,i) (38) 

+ [||u - iSAr|U»(_i,i) + ||u - u»jv||h »(- 1,1)1- 

The interest of this decomposition is that we can now choose tl»j v, still close to u but such that the higher 
norms (i.e. 1)) is uniformly bounded. This will be the case if we choose <b N as being the projection 

of u with respect to the If "(-l, l)-norm over Xn- It is proven in [17] that this element satisfies 

||U - w n \\hH-i,D < cm- ||u|| jrs( -i,i ) , V 7 < o. (39) 
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^From the definition of a and am we have 


(40) 

(41) 


a(xb N , v N ) - a M (viN, «w) = (pd>Ar» , v Na ) - (pdWa, «*»)« 

+ (qw N ,v N )-{qws,v N )u. 

Using Lemma 2.3, first with ip = jnb^ x and then with ip — qw^, we derive that 

|<*(u>Ar,VAr) - <*M(v>N) v fir)\ < C[(2M - 

+ (2 M - ^)--n(M,«)|| 9 ^|,^ i . (M , )( _ M) )| K || ffl( _ M) 

< C[(2M - ^)-™nK-i)||p U ,||^.,., )( _ i i) 

+ (2 M - *)--«(*.') |HI^- (< ... )( _ M)] ||vjv|k l( - M) 

+ C[(2M - N)- min ^-^\\p(u - 
+ (2 M - Ny min ^\\q(u - *i/)|| H ~«(M.o ( _ l(l) )||t»Ar|U.(-x,x) 

Let us recall that min(v, p,a — 1) > 1/2 so that 1) and 1) are algebras. This 

proves that 

|o(tZ»AT,Viv)-aM(tL J v,Viv)| < C((2M- ^)- m<n ^'- 1 )||pu,|| H;?i .(...-0 ( _ lil) 

+ (2 M - iV)- min( '‘-' ) ||g«||^.(,.o ( _ 1>1) ]|Klk.(-i 1 i) 

+ C[(2M - tf)- m< " ( ^- 1 >||p|U ;a .-x,»)||(u - u>xv).||^.c-o ( _ M) 

+ (2 M - N)~ min ^ ||«|U- ( _i.x)ll(« - ^)I| h -»(,.0 ( _ 1i1) ]|K|U‘(-x,x) 

which together with (39) and the fact that M > N (so that 2 M — N > J\T) gives 

KtiW.i,*) - aj,(*#r,»*)| < C[(2M - ^)-'” < "(-'-"- 1 )||pu a || // „ i . ( . > ,-„ ( _ lil) 

+ (2 M - 

+ [iV 1 - ff ||p||tf^(_i,x)||u||Ki(-i,x) 

+ ^'‘ < ’ll9ll«'(-x,i)ll u llw'(-x,x)]|| v iv||w , (-i,i) 

Theorem 2.4 now follows from the fact that if r (—l, 1) is an algebra when r > 1/2. 

Corollary 2,5. Let us suppose that the hypothesis of theorem 2 A are fuMlled and in addition let us suppose 
that M = N .Then the following error estimate holds 

ll« - «jr|Ur‘(-».i) <C[^ m<n <^^ ,0 - l ^>l||p|U i( _i t i ) ||u|| H - j .(, +l ., )( _ 1(X) 

+ Il9|ltf£(-M)NI + II/IIhj ( -x.x)] 

Remark 2.2 — Interpretation of the results: In the practical computations, the more common situation 
is the one where the regularity of the solution is limited by the regularity of the non-constant data. More 
precisely, and especially if we consider the extension to the 2-D case that will be treated in the forthcoming 
section, we shall have o — 1 < min(*/,p, p). In this case, the estimates given in corollary 2.5 prove that 
the optimality of the scheme is achieved with no use of overintegration. This statement is confirmed by the 
numerical experiments shown in Figure 2.1. In the case where the regularity of the solution is better than 
the regularity of the non-constant data, we expect, from theorem 2.4 that some amount of overintegration 
is necessary to recover the optimality of the discretization. The numerical experiments of Figure 2.1 also 
confirm this statement with exactly the rate theoretically proven. It is important to note that the possibility 
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to include some amount of overintegration only is true in the case where the variational formulation is 
the starting point of the discretization. The pure collocation method is incompatible with the concept of 
overintegration. 

2.2 — The Two Dimensional Case. 

Let us suppose here that the problem is set on a square O = (— 1, l) 3 and that no domain decomposition 
is used. The extension of our results to the spectral element case is very simple and it is only for sake of 
limitation of the notations that we restrict our analysis to this particular case. Let us consider now the 
following test problem: Find u E such that 

V . A . Vu = / (42) 

where / is a given force (supposed to be in L 2 ( fi)), A = (a»j) is a 2 x 2 matrix with non-constant coefficients 
and we assume, for sake of simplicity, that it is bounded, symmetric, positive and non-degenerate, i.e. such 
that there exists 2 real numbers p* > 0 and pa > 0 such that 

av • v > Pl (v? + v 3 ), w = (v u v 2 )em 3 (43) 

k,l < Pa i Vi, j = 1,2. (44) 

The variational formulation of this problem : Find u E /fo(O) flUc ^ 

«(*.•)=■((/,»))• Vv E HoW (45) 

will be used as a starting point to define the discrete problem. Here 

= J ' v(0,v>) e (L 3 (n)) 3 , 

v ) = ((A Vu, Vv)), V(u, v) € (H 1 ( n)) a . 

The space of discretization, Xu> consists of all polynomials of degree < N in each variable (ZP//(f])) that 
vanish at the boundary 

x N = Jp N (n)nHl(n) 

The quadrature formula is based on ( M -f l) 2 points derived by a tensor product of the one dimensional 
Gauss-Lobatto formula, and consists of the evaluation Yol x Yql 2 w ^ere the index 1 and 2 refer to the 
first and second spatial direction respectively. The discrete problem is then: Find uu E Xu such that 

aAf(^7V,vjv) = ((/»vat))m, e Xu (46) 

where we have introduced the following discrete scalar product over 2P//(D) 

((*,*))* = £ V(^,^)G(2P W ( fi)) 2 (47) 

QL.XQLa 


and 

aAf(u^,VAr) = ((AVuat, VtiAr))M- 

It is standard to note that, from (43) and (44), this problem has a unique solution and that in the case where 
M = N it is equivalent to a collocation problem baaed on the (N -f l) 3 Gauss-Lobatto points. The analysis 
of the error will be based on Lemma 2.3 and on the following complement of this Lemma 
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Lemma 2.6. Let p be a real number > 1/2. For any V> G //£( — 1, 1), we have 


sup 

v«€V*r 


(V’.vw) - (V’.t'Jv)M 


i’(-l.l) 


< C(2M -N - 1)-" ||^|/Jir»(/0- 


( 48 ) 


* = 1 


Proof: Thifl proof is derived from (29) by taking $ 2 M-n-i instead of $ 2 M-n> which has the effect to cancel 
the middle difference in the right hand side of (29). 


Remark 2.3: It is an easy matter to note that (48) cannot be improved to yield a bound by the best fit by 
polynomials of degree 2M — N as it was in Lemma 2.3. 


We are now in position to give the following error estimate 

Theorem 2.7. Let us suppose that the solution u to problem (42) belongs to H a (Cl), that the coefficients 
(dij ) of the matrix A belong to H^(Q) and that f belongs to H p ( Q) and iet us suppose in addition that the 
real numbers p,cr — l,p are larger than 1. Then the following error estimate holds 

2 

II* - «#f||jr‘(o) < C’[Ar 1 -°||u|| w . ( n ) + (2M - N — £ ||a, y u|| ff » i . ( ^-o ( n) .... 

+ (2M-JSr)^||/||^ (n) ]. 

Proof: The proof is derived from Lemma 2.2 (with obvious changes in the notations) following exactly the 
same lines as in the previous subsection. The only difference relies on the fact that Lemma 2.6 has to be 
used in the analysis of the quantity 


a(u,Vtf) - a M (wM,vx) 

sup — : : : 

||vAr||tf»(n) 


(50) 


Here is a polynomial of degree N (and not N — 1) suitably chosen to be an optimal approximation of u 
in the norm H°{ Q), say, and in any other lower norm from the bidimensional equivalent to (39) (see [16]). 
As in (38) we can write 


sup 


a(tti vn ) - Qjtf (uWi ujv) 


< sup 


a(uj//. Vff) - QAf (u*Af, Vat) 
ll«w|U‘(n) 

+ C||u- lOArlltf.(n). 


The first term on the left-hand side is a sum of four terms that are similar. The first one reads a follows 

((a uw Nx ,v Nx )) - ((auW:rx,v/r x ))Ar 
and can be bounded from Lemma 2.6 as 


((oh^at*, vjvx)) - ((an^Arx.ViVx))^ < C(2M - N - 1) ^||aiiu*|| H mi*( M .r-i)( n )||i;*||x / a(n) 

The reason why (20) is not sufficient to analyze this term is that in the second direction, v^x is of degree 

< N. 

It follows easily from the analysis in the one dimensional case that we can prove 

8up < C (2M - JO-'ll/IUno), (51) 

which allows to derive Theorem 2.7. 
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Corollary 2.8. Let ua suppose that the hypothesis of theorem 2.7 are fullfiUed and let us suppose in addition 
that M = N. Then the following error estimate holds 

||u - ^11^1(0) < C'JV” min ^ l/|ff " 1,p ^{||A||^^(o)||uj|^r(n) + ||/||//'(n)}- 

Remark 2.4 — Interpretation of the results: Ab pointed out in the previous subsection, there is no need to 
overintegrate in the general case where the solution has the same regularity as the non-constant coefficients. 
In 2 dimensions this conclusion still holds, however more weakly. Indeed, we have only been able to prove 
that the plain spectral or spectral element method provides the same accuracy as the best fit by polynomials 
of degree one less. An overintegration with just one more point (M=N-f 1) is sufficient to recover the optimal 
accuracy. Although the results are asymptotically the same, this remark could suggest the use an overin- 
tegration. However, the associated increase in work makes it preferable to use a standard spectral element 
method with polynomial degree of N+l (instead of N) since this is less expensive than an overintegration 
with M=N~f 1 (see appendix). It is only in the case where the spectral method would require much higher 
values of N that overintegration, coupled with the preservation of the tensor product formula, could be of 
interest. Finally, note that the numerical experiments again confirm the theoretical statements with an even 
less important difference between the cases M = N + 1 and M = N. 
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3. ANALYSIS OF SPECTRAL APPROXIMATION IN DEFORMED GEOMETRIES. 

Our purpose is now to extend the previous analysis in the case where the non-constant coefficients 
are induced from the treatment of a non-rectangular geometry. Since the point here is not to analyze the 
behaviour of the approximation when original non-constant coefficients are present in the equation, we shall 
restrict our analysis to the simple problem of a Poisson equation in a deformed domain Q which reads as 
follows: Find u € Hq(Q) such that 

JJ VuVvdxdy = JJ fvdzdy Vv 6 #o(0). (52) 

In order to explain an important issue related to the error analysis, let us first consider the simple case where 
the geometry is rectilinear but not rectangular. 

3.1 — Rectilinear Geometries. 


Let us suppose that the domain O is the trapezoid with vertices A=(0,0), B=(4,0), C=(O t l), and D=(4- 
4/a f l), where a is a real number > 1. It is an easy matter to check that the transformation F : (r, s) h-* (z, y) 
with 


* = K 4 ~ a*)’ 

y- a, 


(53) 


is a one to one maping from the square ft = (0, l) 3 onto the domain 0. We can now define a natural 
correspondance between functions <f> defined over O and functions <j> defined over Q as follows 

4>(r t s) = <j>(F(r t s)). 


This will allow to set the discrete problem. Indeed as we have indicated in the introduction, the spectral 
metod has to be used on square domains. The other domains have to be mapped on such a reference 
square in order to allow for a spectral computation (the spectral element method allows to deal with more 
complex geometries, however, the elements must still be mapped to such a reference square). The variational 
formulation transferred to the domain (l is the following: Find u 6 Hq(Ci) suc ^ 


fvJ(r , a)drds t Vv 6 Hq (ft). 

(54) 

The coefficients J and i = 1,2, 3, are respectively the Jacobian of the transformation F and some smooth 
geometric factor given by 


ff 1 dudv du 


dv dudv , . dudv u ff 

Ts + Ta^ r ' a) ^TrTr W ^ l)]drda = £ 


J = l(a-s), Wl = (l + ^), W3 = (4 — “)(~)i w 3 = (4-^) J . 


(55) 


The discretization of the problem can now proceed and consists in the following: Find its € 2Pw(0)n JSf^(n) 
such th at 

osf(uw.VAr) = ((<7/. ®w))jtfi W € iPw(O) (56) 

where ((., .))u has been defined in (47) and 


,,widu N 

aM(UN,v N ) = {{-j~^r 


dVN 

dr 


))« + (( 


wj dun 

17 


))m 


J dr 

, //Wjduw dv N ^ s , (1)3 du N diff ^ 

+ (( y 17- 17 ))m + (( T 17’ 17 ))m - 


The coefficients in the right-hand side of (56) are extremely regular. We can deduce from the previous 
analysis that no overintegration is required to compute accurately this term in order to achieve optimal 
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approximation. On the contrary, the coefficients that arise in a at on the left-hand side, although analytic, 
may give rise to problem since the factor w*/J has a pole where J vanishes, i.e. at s = a. This pole is not 
in ft but can be very close to it. The similarity of this problem with the one treated in the second section 
(Figure 2.2) seems to imply that, when a is sufficently close to 1, an overintegration is required to compute 
correctly a*/ in order to get an optimal aproximation. However, surprisingly, the experiments performed in 
Figure 3.1 shows that no overintegration is required, even with those values of a that in the one-dimensional 
case induced a big error associated with the consistency error (50). The interpretation of this fact can be 
easily understood, a posteriori, as the pole of 1/J is compensated by the fact that the other part of the 
integrand vanishes and makeB the term 


U)\ flu U>2 du 

Ta7 + Ta7 + '" 

more regular than it appears. In order to be more precise, let us state 

Lemma 3.1. Let us assume that the solution u to problem (52) belongs to H° (£1) for some real number 
a > 2. Then we have 


sup 

0 n€JPn(Ci) 


°(fl| ^n) - <*M(^Ar»%) 
II*MIh»(A) 


< C[||u - u>w||„ l(A) + N l ||4 - 

+ (2M -N -l) 1 


, du, 


+ 11^1 


(57) 


a x "«'-*(A)' r|l ( 9y ll "'- , (A) 


J* 


Proof: Let us go back to the original equation (52) and remark that 


L 


VuWvdxdy 




du w 4 »\ , &vn , w „ 4 

-£7)( r ,') + ^( r *')[(-57( r >')(-) + -5T( r «')( 4 " -•Wrds 


and that a simple change of variable formula yields 

(( 


(((4 ~ “*) 


flu 

dv N 

d x' 

dr 

du 

dv N 

V 

dr 

du 

dv N 

V 

ds 


xx _//!£* a *H \\ 

"J dr' dr 


nt w \™ ov N S s _ut u dfl Afxx , ,<*2°* ov N 

((( T } %' ~97 ))m - ~a^ ))M + (( 7al' "sr ))M 


b>2 0U dVjf , 


w . , ,,<*>$ db flthv xv 

))m — (( j ))m + (( 7 » Rm ))M 


J da' da 


so that 
where 


a(u, v N ) — dM (wn, vn) = A\ + Bi + A2 + £2,1 + £2,2 + A3 + #3,1 + #3,2 


Ax 




d ” N \t \A A ff dU 

—)(r,s)drds ))„ 


Bl ~ “l~dr' ~d7 ))u ~ (( J IT* ~d7 ))M 
a ff au f , u 4r.. ... 4r.5u 

M = lit, ~ «<—>*• '»r )) " 

d -m w di>N u ,, ui-l dw N dv N 
B it i - (((— 7— )-5r, — ))m - (((— 7—)-^-. 


*2,2 = (((=£) 


dr' dr 
U)2 x du dVN 


J } da* dr 


xx ((( U2.dw N dv N 


(58) 
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M = Jfa ^ (r ' ,)( Tr (r ' ,)(4 - i*)** - (( < 4 - TT»" 

o ~{(( U 3\ d{l diN \\ <(( U3 \ d * N diiN \\ 

B ** - > t*' ~ar ))u ~ ((( y )_ 5T’ yr }) "- 

In the decomposition (58), the terms A _ are small from the accuracy of the quadrature formula; the terms 
2?_ are small since ws is a polynomial that is close to u. More precisely, let us denote by C the supremum 
of \l/J\, 1(1 - Wi)/J\ t \u) 2 /j\ 1 \u> 3 /J\ } on H. It is an easy matter to derive that, for any index a, where 
a = 1, (2, 1), (2, 2), (3, 1), (3, 2) we have 

\B a \ < C\\I M ^{u - ^/v)|| L a(A)||^|| // i(A), 
where t in the derivation stands for a or r. Using (33), with p = 1, we derive 

|#a| £ ^[||6“ u^II//*(A) + N 

< CfHu - + N 1 ||ti - llj/i(A)* ( 5 ®) 

Let us denote now by c the supremum of 1, |(4 — 4s/a)|, |4r/a|, on ft. We derive as in section 2.2 that 

\Ai\ < c(2M - N - 1) 1 ~ <7 [||^Hh^-»(A) + ll^lli/'- l (A)]ll^llir l (A) 

and the Lemma follows from (58), (59) and (60). 

In order to analyse the error in the approximation of problem (52) with (56), we first note that (16) 
(with obvious changes of notations) also holds in this case. This leads to the following: 

Theorem 3.2. Let us suppose that the solution u of problem (52) belongs to H° ( O) and that / belongs to 
H p ( ft) and let us suppose in addition that the real numbers a - 1 ,p are larger than 1. Then the following 
inequality holds 

||U - U* ||„, (A) < C[W l - a |Mltf' ( n ) + C(2M -N- l) 1 -||u||*.<n) 

+ C( 2M - N)-'\\f\\ H , {(1) 

Proof: This is an easy consequence of (48) and of Lemma 3.1, if we choose wn as being the best fit of u in 
H a (Cl) by elements of 1P^(Q). 


(60) 


(61) 


3.2 — General deformations. 

Let us suppose in this subsection that the domain ft is the image of fl by some mapping F : (r, s) h-* 
(x = Fi(r, s),y = F 2 (r, s)). We no longer assume that F is bilinear, nevertheless we assume it has some 
regularity, more precisely, that it is in H^(Q) for some p > l, and that its Jacobian J is larger than some 
Pi > 0. The problem (52) is transferred onto H as before and reads : 

Find u G Hq(Q) such that (54) is satisfied, 
where the geometric factors i = 1,2, 3 are defined as follows 

A 

LOi = (- 




da 


r dFi dF\ dF 2 dF 2 , 

“ “l o- a ~ 1 


da dr da dr 
.. _ ( dFl \* + ( dF *\ 2 


( 62 ) 


The discrete problem is now: 

Find un £ 2Pjv(ft) n Hq(Ci) such that (56) is satisfied . 

Using here the same technique as in the proof of Lemma 3.1, and theorem 2.4, we obtain 
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Lemma 3.3. Let us suppose that the solution u to problem (52) belongs to H° { ft), for some real number 
a > 2, that the deformation mapping F belongs to for some real number p > 1. Then the following 

inequality holds 


+ (2M — N — IU. (ft) |||^||J r _ i . ( + llglli-„....„ (0) ]) 


(84) 


Remark 3.2: It is important, at this stage, to note that the regularity of u and the regularity of where 
t is either s or r is bounded by the regularity of both u and F, which means that 




(65) 


is the best result we can expect. Note also (but this U natural after the analysis of section 3.1) that it is 
only the regularity of F that is important and not at all the regularity of its inverse (as one could previously 
fear from the factor 1 /J). 


iFrom Lemma 3.3 and (65), we derive in a (now) standard manner the 

Theorem 3.4. Let us suppose that the solution u belongs to problem (52), belongs to H°(Q), for some real 
number a > 2, that the deformation mapping F belongs to // M (Q) for some real number p > 1 and that the 
forcing function f belongs to H P (Q) for some real number p > 1. Then the following inequality holds 

lift - < C{Ar- m -('‘-^- ‘)||F||^ (A) ||u||^. (< .. 0(0) 

+ (2 M-N- (66) 

+ (2 M - /)(n)} 

Remark 3.3: In the light of the previous theorem, it appears that there is no need to overintegrate in 
the case where the deformation is responsible for non-constant coefficients in the formulation of the discrete 
problem, even if the geometry is very distorted. This is readily seen by comparing the results obtained in the 
theorem by choosing M » N and M = N (the gain of just one degree is not, as explained in the remark 
2.4 sufficient to start the overintegration machinery). We can explain it a posteriori now as overintegration 
could allow to improve part of the consistency error term (57) in the case where the geometry mapping F is 
less regular then the solution u itself ( contrary to the case in section 2, this can happen, especially if domain 
decomposition is used). However, as pointed out in the previous remark, loss of regularity of F induces also 
a loss of regularity of ti which is the function that is of importance since it is this one that is approximated 
by polynomials. 


Remark 3.4: We have not considered here the case where the geometry mapping F itself is approximated 
by isoparametric polynomials. This effect can easily be analysed by following the same lines as previously 
and including the standard arguments adapted to the finite element discretization [18] 
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APPENDIX (in collaboration with A. T. PATERA.) 

In this appendix we consider the computational cost of using overintegration. We show that the matrix- 
vector products can still be evaluated efficiently using tensor product forms, and that for a fixed discretization 
(polynomial degree N) t the increase in operation count is at least a factor of d in IR d . Overintegration is 
attractive to use only if the convergence rate (error as a function of polynomial degree) improves substantially. 
In this case we can achieve a fixed accuracy either with a low polynomial degree and overintegration, or a 
higher polynomial degree with consistent quadrature. 

With no loss of generality we consider here the evaluation of the left hand side of (46), 

<Ur (ujv> vn) = ((pVujv Vvjv))a/ (67) 

Here p represents a non-constant coefficient which either comes from the original strong formulation (as in 
(42)), or represents a geometric factor in the case of deformed geometry, or both. The evaluation of (67) 
represents the computationally most expensive part in an iterative solution of (46). 

In order to implement (46) we require a basis for our high-order polynomial space Xs. The choice 
of basis does not effect the error estimates, however it greatly effects the conditioning and sparsity of the 
resulting set of algebraic equations, and is critical for the efficiency of parallel iterative solution procedures. 
We choose an interpolant basis to represent w n E 

N N 

( 68 )> 

p—0 q~0 

where (z,y) E 0 — ► (r,s) €] - 1, 1[ 2 . Here the h p (z) are the one-dimensional N th - order Lagrangian inter- 
polants through the Gauss-Lobatto Legendre points ( p (hp E IFV(] — 1 , 1 [)» (^ ) — and w P q is the 
value of ws at the local node {(p,£q)- In addition to (68) we require w n E Xn to honor the homogeneous 
Dirichlet boundary conditions. 

The bases (68) are now inserted into (67), and the test functions are systematically chosen to be unity 
at one global node and zero at all the other Gauss-Lobatto Legendre points. We then arrive at the following 
discrete statement, 

AT Af N N 

a Af i^N) v n) = ^ y ] PaPpPapDailfJj DamlpnVmn + (69) 

a=0 fizzO rn — On — O 

where we for convenience only have included the derivatives with respect to the first (local) spatial direction. 
Here p v are the one-dimensional quadrature weights, and D and / are the one-dimensional derivative and 
interpolation operators defined as: 


dh l± 

Dpt — ^ (fp)i 

(70) 

Q* 

1! 

o* 

(71) 


Using tensor-product sum-factorization techniques [5] we can now evaluate (69) efficiently as follows: 

{D a i(Ipj (papp (Pa/?(^om(f^n u mn))))))i 

where the expression in the innermost partenthesis is evaluated first. An evaluation in this order results in 
the following operation count, 

MN 2 + M 2 N + M 2 + M 2 + M 2 N + MN 2 = 2(M 2 N + MN 2 + Af 2 ). 

In the case where M = N (no overintegration), the interpolation operator becomes J pq = 6 pq and (69) 
reduces to the following expression, 

M N 

» V/v) = y ^ y^ PaPjPaj D a iD arn Urnj H" •••• 

a=0 m=0 


ifl 


which can be efficiently evaluated in 2 (N 3 + N 2 ) operations. 

We are now in a position to compare the computational cost with and without overintegration. In the 
case M = JV + 1 the operation count is approximately twice the cost of using consistent quadrature (Af = N ) 
(in three dimensions the cost would increase with a factor of three). If M = 2 N } say, the cost to evaluate 
(67) is four times the cost with M — N. The factors which determine whether overintegration is economical 
are the the difference in convergence rate with and without overintegration, and the specified accuracy of 
the discrete solution. 
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FIGURES 


Figure 2.1. A plot of various numerical results relevant to the approximation of the one dimensional problem 
(1), in order to analyze the effect of overint eg ration. We have considered here the case where q = 0 and 
different values of p and / are simulated and the effect of overintegration is adressed. The relation between 
p and / are always chosen in such a way that the exact solution is u = sin(3xx)e~ s (of course, we do not 
use the knowledge of the exact solution in the simulations 1!). The discretization is based on if = 4 equal 
spectral elements for different values of the polynomial degree, N. First of all, since this is the aim of the 
numerical scheme, we have computed the best fit of u in Hq(Q> 1) by piecewise polynomials of degree < N 
(plot ▲); the plot proves a convergence better then exponential that is consistant with the fact that the 
solution is entire. Next we consider the case p = 1 with no overintegration (plot o) and we can check that 
the approximation is optimal and very close to the best fit. 

We then choose p = for various values of a. The case where a = .25, treated with no overintegration, 
is represented on the plot 0 and we can remark that the approximation is still very good. The convergence 
rate does not however follow the best fit of u, but the best fit of p (consistently with the corollary 2.5) that is 
a straight line in this scale since p has a pole outside (0,1). Let us now use a more singular p, corresponding 
to a = .05. The plot 0 represents the results when no overintegration in used. It is very close to three other 
plots representing the results when overintegration (Af = SN/2) is used to compute only the right-hand 
side (the contribution of /) (□) or only the left-hand side (corresponding to p) (plot A), while the third 
one represents the best fit of p (plot •) in the L 2 ( 0, l)-norm by the discrete functions. In order to recover 
optimal results, overintegration must be used on both the left and the right-hand side of the equation (plot 
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ERROR 


Figure 2.2. This Figure represents the equivalent of the previous case in a two dimensional domain. The 
equation is here (42) with A = pld and different values of p and /, computed so that the solution is always 
u = jin(ffx) y(l — The domain of computation fi is (0, l) a and decomposed in four squares. The plot 

• represents the case p = 1 with A/ = Af, A/ = AT + 1, M = 3 N (there is no difference between the two last 
experiments and the first case is a bit less precise but cannot really be distinguished from the two last cases). 
The plot o represents the experiment with p = that is a very smooth function and a treatment with 
no overintegration gives the same accuracy as the previous plot (that certainly also corresponds to the beat 
fit). Now comes the treatment of the case p = 1 + cos{2nx + 3xy) which is a smooth function but is worse 
than u. The plot ▲ corresponds to no overintegration and we clealy see that some degree of overintegration 
is necessary to recover the optimal accuracy, as shown in the case where M = Af + 1 (plot □) and M — 3 N 
(plot A). 

Figure 3.1. On this plot we prove that the singularity of the function i* not relevant for the need 

of overintegration. We consider the problem is (52) for two different domain decompositions A and B. They 
are 


A 



The three plots represent the discretisation error in the semi-norm with no overintegration. The plot o 
corresponds to the solution u = *in(2y) for the two cases A and B. They both coincide since the deformation 
is in the x-direction and the solution only depends on y. Next, we repeat the experiments with u = stn(|) 
and we arrive to the (a priori) surprising result: plot A for A, and □ for B. We refer to section 3.1 for the 
explanation. The fact that B is even better in this case is due to the fact that the discretisation on the two 
domains is more equilibrated and that the left domain in case A is too long. 



Figure 2.2 
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Figure 3.2. Thu plot is to prove that the overintegration can even have a negative effect. Again we consider 
equation (52) for p = 1, with a decomposition of n = [0,7] x [0, l] into two subdomains: 



The interface is described by the equation x = j - ||y| 3 . This equation is not regular and its interpolant 
is more rough when the degree is higher. Since this interpolant is the only geometric factor for the discrete 
problem we can understand the plots A and ▲ that correspond to M =* N and A/ = 2 N t respectively, for 
the approximation of the solution u = *tn(|). As before, if the deformation and the solution are in different 
directions (this is a very particular case 111), for example u = stn2y, the behaviour of the approximation is 
particular also and this is the only case where overintegration can help (in the case of deformed geometries) 
as is illustrated in the plots o and • corresponding to M = N and M = 2 N respectively. 

Figure 3.3. A plot of the discretisation error in the semi-norm when solving (52) with p = 0. The domain 
fl = [0, 2] x [0, l] is decomposed into two subdomains: 


and this time the geometry is tortured in the following way. The subdomain on the left is mapped onto the 
reference square by the mapping x =* r + 0.5(|r| — l)r 3 , while the mapping is affine in the other subdomain. 
The plot A correspond to a solution u = *in(|) that only depends on x and is treated both with M = N 
and M = 2 N. The plot o is corresponds again to a particular case as it is related to u — *tn(2y) in both 
cases M — N and M = 2N. 



Figure 3.2 
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